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Abstract 

Building on substantial foundational progress in understanding the 
effect of a small body’s self-field on its own motion, the past 15 years 
has seen the emergence of several strategies for explicitly computing 
self-field corrections to the equations of motion of a small, point-like 
charge. These approaches broadly fall into three categories: (i) mode- 
sum regularization, (ii) effective source approaches and (iii) worldline 
convolution methods. This paper reviews the various approaches and 
gives details of how each one is implemented in practice, highlighting 
some of the key features in each case. 


1 Introduction 

Compact-object binaries are amongst the most compelling sources of grav¬ 
itational waves. In particular, the ubiquity of supernrassive black holes 
residing in galactic centres jjf has made the extreme mass ratio regime a 
prime target for the eLISA mission [0-0] • Meanwhile, the comparable and 
intermediate mass ratio regimes are an intriguing target for study by the 
imminent Advanced LIGO detector [fj]. In order to maximise the scientific 
gain realised from gravitational-wave observations, highly accurate models 
of gravitational-wave sources are essential. For the case of extreme mass 
ratio inspirals (EMRIs) — binary systems in which a compact, solar mass 
object inspirals into an approximately million solar mass black hole — the 
demands of gravitational wave astronomy are particularly stringent; the 
promise of groundbreaking scientific advances — including precision tests 
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of general relativity in the strong-field regime @-03 and a better census of 
black hole populations - hinges on our ability to track the phase of their 
gravitational waveforms throughout the long inspiral, with an accuracy of 
better than 1 part in 10,000 [2]. This, in turn requires highly accurate, long 
time models of the orbital motion. 

For the past two decades, these demands have stimulated an intense pe¬ 
riod of EMRI research among the gravitational physics community. Despite 
of the impressive progress made by numerical relativists towards tackling 
the two-body problem in general relativity (for reviews see Refs. HM1), 
the disparity of length scales characterising the EMRI regime is a signifi¬ 
cant roadblock for existing numerical relativity techniques. Indeed, to this 
day EMRIs remain intractable by current numerical relativity methods, and 
successful approaches have instead tackled the problem perturbatively or 
through post-Newtonian approximations (see Ref. 14] for a review). This 
article will focus on the first of these; by treating the smaller object as a 
perturbation to the larger mass, the so-called “self-force approach” reviewed 
here has been a resoundingly successful tool for EMRI research. 

Within the self-force approach, the smaller mass, fi « M, is assumed 
to be sufficiently small that it may be used as a perturbative expansion pa¬ 
rameter in the background of the larger mass, M. Expanding the Einstein 
equation in /r, we see the smaller object as an effective point particle gen¬ 
erating a perturbation about the background of the larger mass. At zeroth 
order in /r, the smaller object merely follows a geodesic of the background. 
At first order in /r, it deviates from this geodesic due to its interaction with 
its self-field. Viewing this deviation as a force acting on the smaller ob¬ 
ject, the calculation of this self-force is critical to the accurate modelling 
of the evolution of the system. For the purposes of producing an accurate 
waveform for space-based detectors, and for producing accurate intermedi¬ 
ate mass ratio inspiral (IMRI) models, it will be necessary to include further 


effects up to second perturbative order [15], [16|. Indeed, recent compelling 
work 17, suggests that IMRIs - and even comparable mass binaries - 
may be modelled using self-force techniques. 

A naive calculation of the first order perturbation due to a point particle 
leads to a retarded field which diverges at the location of the particle. The 
self-force, being the derivative of the field, also diverges at the location of 
the particle and one obtains equations of motion which are not well-defined 
and must be regularized. A series of formal derivations of the regularized 
first order equations of motion (now commonly referred to as the MiSa- 
TaQuWa equations, named after Mino, Sasaki, Tanaka [191 ] and Quinn and 
Wald 20] who first derived them) for a point particle in curved spacetime 
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culminating in a rigorous work by Gralla and 
in the gravitational case and by Gralla, et al. 
in the electromagnetic case. This was subsequently extended to second 
perturbative order by Rosenthal 34j-f37], Pound 38 h 40|], Gralla 41] and De- 
tweiler 42]. These derivations eliminate the ambiguities associated with the 
divergent self-field of a point particle and provide a well-defined, finite equa¬ 
tion of motion. Building upon this foundational progress, several practical 
computational strategies have emerged from these formal derivations: 


• Dissipative self-force approaches: While the full first-order self-force 
is divergent, it turns out that the dissipative component is finite and 
requires no regularization. This fact has prompted the development of 
methods for computing the dissipative component alone, sidestepping 
the issue of regularization altogether. These dissipative approaches 
fall into two categories: 


1. 


Flux methods: By measuring the orbit-averaged flux of gravita¬ 
tional waves onto the horizon of the larger black hole and out 
to infinity, the fact that the field is evaluated far away from the 
worldline means that no divergent quantities are ever encoun¬ 
tered. This approach yields the time- average^ dissipative com¬ 
ponent of the self-force 4^- 5 


2. Local/instantaneous dissipative self-force: The time averaging el¬ 
ement of flux methods can be eliminated by instead computing 
the local, instantaneous dissipative component of the self-force 
from the half-advanced-minus-half-retarded field 43, 51-531]. 


Both methods, however, fundamentally rely on neglecting potentially 
important conservative effects which can significantly alter the orbital 
phase of the system. 


The mode-sum approach: Introduced in Refs. |54|, 550, and having 
since been successfully used in many applications, the approach re¬ 
lies on the decomposition of the retarded field into spherical har¬ 
monic modes (which are finite, but not differentiable at the parti¬ 
cle), numerically solving for each mode independently and subtracting 
analytically-derived “regularization parameters”, then summing over 
modes. 


1 For the case of inclined orbits in Kerr spacetime, this is more appropriately formulated 
as a torus-average. 
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The effective source approach: Proposed in [56| and 57], the approach 
implements the regularization before solving the wave equation. This 
has the advantage that all quantities are finite throughout the calcu¬ 
lation and one can directly solve a wave equation for the regularized 
field. 


The worldline convolution approach: First suggested in [58j, |oyj, one 
computes the regularized retarded field as a convolution of the retarded 
Green function along the past worldline of the particle. Although the 
approach is the most closely related to the early formal derivations, it 
is only recently that it has been successfully applied to calculations in 
black hole spacetimes. 


For a comprehensive review of the self-force problem, see Refs. }6Cn63l] . In 
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this paper, I will review the various approaches and give details of how each 
one is implemented in practice, highlighting the advantages and disadvan¬ 
tages in each case. 

This paper follows the conventions of Misner, Thorne and Wheeler 
a “mostly positive” metric signature, (—, +, +, +), is used for the spacetime 
metric, the connection coefficients are defined by = \g Xa + 9 av^ ~ 
dnv.a), the Riemann tensor is R a Xflu = - T^T^, 

the Ricci tensor and scalar are R a /3 = R ^ a ^3 and R = R a a , and the Ein¬ 
stein equations are G a p = R a p — \g a gR = 87t T a p. Standard geometrized 
units are used, with c = G = 1. Greek indices are used for four-dimensional 
spacetime components, symmetrisation of indices is denoted using parenthe¬ 
sis [e.g. (cc/3)], anti-symmetrisation is denoted using square brackets (e.g. 
[a/3]) and indices are excluded from symmetrisation by surrounding them 
by vertical bars [e.g. (a|/3|7)]. Latin letters starting from i are used for 
indices summed only over spatial dimensions and capital letters are used to 
denote the spinorial/tensorial indices appropriate to the field being consid¬ 
ered. Either x or are used when referring to a spacetime field point and 
z(t ) or z^(t) are used when referring to a point on a worldline parametrised 
by proper time r. Finally, a retarded (or source) point is denoted using a 
prime, i.e. z'. 


2 Equations of Motion 

The formal equations of motion of a compact object moving in a curved 
spacetime are now well established up to second perturbative order. Writ¬ 
ing the perturbed spacetime in terms of a background plus perturbation, 
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g a p = g'a'p + h a p, the equations of motion essentially amount to those of 
an accelerated worldline in the background spacetime, with the acceleration 
given by a well-defined regular field which is sourced by the worldline. To 
order g, this coupled system of equations for the worldline and its self-held 
are commonly referred to as the MiSaTaQuWa equations and are given (in 
Lorenz gauge, assuming a Ricci-flat background spacetime) by 

+ 2C a y / 3 5 h™s = —167r/r J g a >( a u a ' gp)pvP'8±{x,z(T'))dT' (la) 

ga a = gk a ^ s h^. 5 (lb) 


with 


k a^8 = 1- g^u'u 6 - iu'VuV + \u a g^ )U s + \g^y 


Here, g is mass of the object, g a ’ a is the bivector of parallel transport, 
C a / 3 is the Weyl tensor of the background spacetime, and we use the trace- 
reversed metric perturbation h a p = h a p — h = h 7 7 . 

One can also consider compact objects possessing other types of charge. 
For example, a particularly simple case is that of a scalar charge q with mass 
m and scalar held <f>, in which case the equations of motion are given b}HH 


(□ - £R)<f> ret 


dm 

dr 


—47t<7 J 64(x,z(r'))dT' 

(2a) 

q(g^ ) +u a u^)^ 

(2b) 

-qu a <S> R a . 

(2c) 


Here, R is the Ricci scalar of the background spacetime and £ is the coupling 
to scalar curvature. Similarly, for an electric charge, e, one obtains equations 
of motion which are given in Lorenz gauge by 


□H^ et - RaAJ 1 = -47re J g aa 'U a ’5 4 (x,z(t')) dr' (3a) 

ma a = e(g^ ) + v^u^A^^u 1 , (3b) 

2 In the scalar case, it is important to distinguish between the self-force F a = V Q< f> and 
the self-acceleration, which is given by projecting to self-force orthogonal to the worldline, 
o“ = {9(0) -\-u ol u^)Fp. In the electromagnetic and gravitational cases the self-force has no 
component along the worldline and the two may be used interchangeably. 

3 We assume that the mass m is small and ignore its effect on the equations of motion. 
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where R a p is the Ricci tensor of the background spacetime and A^ is the 
vector potential. 

The key component in all instances is the identification of the appropri¬ 
ate regularized field on the worldline. Detweiler and Whiting identified a 
particularly elegant choice for the regularized field, written in terms of the 
difference between the retarded field and a locally-defined singular field, 

$ R = $ ret_ $S) A R = A ret _ ^ h^, = h$ - h^. (4) 

In addition to giving the physically-correct self-force, the Detweiler-Whiting 
regular field has the appealing feature of being a solution of the homoge¬ 
neous field equations in the vicinity of the worldline. Most computational 
strategies essentially amount to differing ways of representing this singular 
fieltfl and obtaining the regularized field on the worldline. 


3 Numerical regularization strategies 

In a numerical implementation, it is essential to avoid the evaluation of 
divergent quantities. In the case of self-force calculations, both the retarded 
and singular fields diverge on the worldline so one must avoid evaluating 
them there. Several strategies for doing so have emerged over the years (see 
Table |T] for a summary). 

One option is to only ever evaluate finite, dissipative quantities (e.g. 
the retarded field far from the worldline or the half-advanced-minus-half- 
retarded field on the worldline). This is the basis of the dissipative methods 
mentioned in the introduction. Since these methods effectively avoid the 
problem of regularization, they will not be discussed further here; we will 
return to them in Sec. [4j It is worth noting, however, that the methods typ¬ 
ically used by dissipative calculations are essentially the same as those used 
by mode-sum regularization for computing the retarded field, but without 
the additional regularization step. 

This leaves three regularization strategies which allow the regularized 
field to be computed on the worldline without encountering numerical di¬ 
vergences: worldline convolution, mode-sum regularization, and the effective 
source approach. 

4 Some methods [§5],[66) rely on alternative prescriptions for the singular field than that 
proposed by Detweiler and Whiting. 
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Case 


Worldline 


Mode-sum 


Effective Source 


c3 

o 

to 


rd 

V 

CD 

N 

c3 


CO 


circular (apprx) 
generic 

(quasilocal) M S; 
generic 169 Tl|]; 
static [t^: 
accelerated 


radial FFSl : 
ircular [75 


circular [75-781; 


eccentric 
static 


72fl; 


79 8 . 


circular 56, 57, 65, 

M|; 

eccentric 18711 : 
evolving [881 ]; 


Sh 
P-1 
CD 

W 


generic [68l |; 
accelerated 


circular [891; 
equatorial 9(3, Ell];. 
inclined circular [92j] ; 
accelerated 93f]; 


circular [961; 


static 


94, 9 


eccentric 


97| ; 


H 


13 

CD 

N 

Sh 

03 


HH 

U 

CO 


static 7: 




static 


eccentric [82|, [98|; 
static (Schwarzschild- 
de Sitter) [99| ; 
radial (Reissner- 
Nordstrom) Hit]: 


f—l 

f—( 
CD 


equatorial [90]; 
accelerated 95 


I 

0 


13 

CA 

N 

Cw 


CD 

CO 


generic 

(quasilocal) [lQl ]: 


radial [] 
circular 
eccentric 


lqiujJi 


wm- 


circular 122 ]; 


12QD; 

osculating 


1211 ]; 


CD 

W 


circular 
(quasilocal) 151 
branch cut 1125 


equatorial 9 
accelerated 


circular 119f, 124]; 


circular 125); 
generic 126i] : 


Table 1: Summary of regularization methods employed by self-force calcu¬ 
lations in black hole spacetimes. 
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3.1 Worldline convolution 


The worldline convolution method relies on a split of the regularized self¬ 
force into an “instantaneous” piece and a history-dependent term. The 
instantaneous piece is easily calculated from local quantities evaluated at 
the particle’s position, 



The history-dependent term is much more difficult to calculate, as it is given 
in terms of a convolution of the derivative of the retarded Green function 
along the worldline’s entire past-history (see Fig. [1]), 


$ h i st = q f V Q G ret [x’, z{t')\ dr ', (6a) 

J — OO 

4$ = e f [x, z(t')]u<*' dr ', (6b) 

J — OO 

= ^ [ ^i G % a 'p'[ x ->z(T')\u a 'uP' dr'. (6c) 

The covariant derivatives here are taken with respect to the first argument 
of the retarded Green function. Likewise, the regularized self-field can be 
obtained from a worldline convolution of the retarded Green function itself; 
similar formulae can also be derived for higher derivatives. The retarded 
Green function appearing in these equations is a solution of the wave equa¬ 
tion with a delta-function source, 

(□ - £R)G ret = -4vr <5 4 (z, x'), (7a) 

□<w - R a ^Gf a , = -4;r W <5 4 (x, x'), (7b) 

+ 27? Q 7 i a <5 Gl ) ,5 t a / ( g/ = — 4:TT g aa igppi5 4 (x, x'), (7c) 

with boundary conditions such that the solutions correspond to purely out¬ 
going radiation at infinity and no radiation emerging from the horizon. 

The regularization (i.e. subtraction of the Detweiler-Whiting singular 
field) is formally achieved by the limiting procedure in the upper limit of 
integration, i.e. by cutting off the integration at t~, slightly before the 
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Timelike worldline of the particle 
Current location of the particle - z(z) 

Integral in quasilocal region 

/ V v G ret {z{T),z(T'))dT' 

Matching point - z(z m ) 


Boundary of causal domain where 
Hadamard form is valid 

First null geodesic which 
re-intersects the worldline 


Integral outside quasilocal region 

I '" X7 a G ret ( Z (T),z(T'))dT' 
J — oo 


Figure 1: Schematic representation of the worldline convolution method. 
The equations shown are for the case of the self-force on a scalar charge, 
but are representative of similar equations in the electromagnetic and grav¬ 


itational cases. Reproduced from Ref. [69 ] 


coincidence point z = x. In practice, this is done by taking the integration all 
the way up to coincidence, but excluding the direct contribution to the Green 
function at coincidence (i.e. the term proportional to 6(a) in the Hadamard 
form for the Green function). Although the retarded Green function also 
diverges at certain other points along the past worldline (in particular at 
null-geodesic intersections where the particle sees null rays emitted from its 
past), it turns out that these are all integrable singularities (of the form 1/a 
and 6(a)) and the integral may be accurately evaluated to give a finite and 
accurate value for the regularized self-force. 

Despite having been proposed in the early days of EMRI self-force cal¬ 
culations 58], the worldline convolution method was largely ignored for a 
long time by numerical implementations (the notable exception being in¬ 


vestigative studies [59|, 122] which did not complete a full calculation of 
the self-force). A likely reason is that the method relies on knowledge of 
the retarded Green function for all points on the particle’s past worldline. 
While methods have existed for decades for computing portions of the re- 
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tarded Green function, it turns out to be very difficult to obtain it accurately 
everywhere that it is needed for the worldline convolution. 

Thankfully, recent progress has led to two practical computational strate¬ 
gies, both of which have been successfully applied to compute the self-force 
in black hole spacetimes. The first of these is based on a frequency-domain 
decomposition of the Green function, and builds on the rich history of black 
hole perturbation theory developed over the past several decades. The sec¬ 
ond, a time-domain approach, has been a very recent development and shows 
a great deal of promise for the future. 

Independently of whether a frequency-domain or a time-domain scheme 
is used, a common problem with both is that they fail at early times when 
source and field points are close together; in the frequency-domain case the 
convergence is poor at early times, while in the time-domain case features 
from the numerical approximation pollute the data at early times. A rela¬ 
tively straightforward solution which has been successfully applied in both 
scenarios is to only rely on their results at late times and to supplement 
them at early times with a quasilocal Taylor series expansion. Provided 
a sufficiently early time can be chosen where both the distant past and 
quasilocal calculations converge, this yields a global approximation for the 
Green function which is sufficient for producing an accurate result for the 
self-force. To date this has been shown to be possible in the case of a scalar 
charge in Nariai 128;], Schwarzschild [fit], 70] and Kerr [ 129 ] spacetimes. 


Given an approximation for the Green function valid throughout the past 
worldline, it is then trivial to numerically integrate Eq. ([6]) to obtain the self¬ 
force (see Fig. [I]). As mentioned previously, the divergences in the retarded 
Green function at null-geodesic intersections on the past worldline are all 
integrable singularities and do not pose a significant obstacle to accurate 
numerical evaluation of the integrals. 


3.1.1 Quasilocal expansion 

In order to obtain an approximation to the retarded Green function which 
is valid at early times, it is convenient to start with the Hadamard form for 
the Green function (in Lorenz gauge) 


G ret 
AB 


( x,x ') = ©_ Uab' 5 ((j) - Vab'Q(-v) 


( 8 ) 


where ©_ is analogous to the Heaviside step-function, being 1 when x' is 
in the causal past of x, and 0 otherwise, 8(a) is the covariant form of the 
Dirac delta function, Uab' and Vab' are symmetric bi-spinors/tensors and 
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are regular for x' —> x. The bi-scalar a ( x , x') is the Synge world function, 
which is equal to one half of the squared geodesic distance between x and 
x'. In particular, a(x,x) = 0 and a(x,z) < 0 when x and z are timelike 
separated. Because of the limiting procedure in the history integral, Eq. (J 6 |) , 
only the term involving Vab' is non-zero, and at all required points along 
the worldline 0_ = 1 = ©(— a). The problem of determining the retarded 
Green function at early times therefore reduces to finding an approximation 
for Vab'(x,x') which is valid for x and x' close together. 

Several methods have been developed for computing approximations to 
Vab 1 ■ Fundamentally, they rely on either the use of a series expansion, or 
on the use of numerically evolved transport equations (ordinary differential 
equations defined along a worldline). The series expansion approach has 
been the most fruitful to date with results including: leading-order coor¬ 
dinate expansions in Schwarzschild and Kerr spacetimes for scalar 


681 , 


and gravitational cases 101 |], high-order coordinate expansions in spheri¬ 
cally symmetric spacetimes (including Schwarzschild) [67(, formal covariant 
expansions in generic spacetimes 130^ 133 ], and moderately high-order coor¬ 
dinate expansions in Schwarzschild~| 8 ^~and Kerr |Qo| spacetimes. The only 
numerical calculation I am aware of was done in [1331 ] for generic spacetimes 
(with an example application in Schwarzschild spacetime). 

The series expansion method produces an expression for V(x,x') as a 
power series in the coordinate distance between x and x'. For example, for 
the scalar case in Schwarzschild spacetime it takes the form 


V{x,x') = Y v ijk(r) {t - - cosj) 3 (r - r') k , 


(9) 


i,j,k=0 


where 7 is the angular separation of the points and the Vijk are analytic 
functions of r and M. It is straightforward to take partial derivatives of 
these expressions at either spacetime point to obtain the derivative of the 
Green function. Although this series on its own may be sufficient for use in 
the quasilocal component of a worldline convolution, it turns out that some 
simple tricks allow for a vast improvement in accuracy. It turns out that, be¬ 
cause V(x,x') diverges at the edge of the normal neighbourhood, the series 
approximation benefits significantly from Pade resummation which incorpo¬ 
rates information about the form of the divergence [67j. One minor caveat 
is that since Pade re-summation is only well defined for series expansions in 
a single variable, it is necessary to first expand r' and 7 in a Taylor series 
in t — t ', using the equations of motion to determine the higher derivatives 
appearing in the series coefficients. Then, with V(x,x') written as a power 
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series in t — t' alone, a standard diagonal Pade approximant provides an 
accurate representation of the Green function in the quasi-local region. 


3.1.2 Frequency domain methods 


Frequency domain methods for computing the retarded Green function rely 
crucially on the separability of the wave equation. In the scalar, Schwarzschild 
case that represents the current state-of-the arid [69] (also see fliil for a re¬ 
lated calculation in Nariai spacetime), this can be achieved by writing the 
Green function as a sum of spherical harmonic and Fourier modes 


-| ^ | poo-\-te 

( X,x') = —7^PHcos7)— / 

rr £=0 27r J-oo+ic 


Qret / ^ _/ 


( 10 ) 


Here, e > 0 is a formal parameter to ensure the correct boundary conditions 
are satisfied for a retarded Green function. Substituting this into the wave 
equation, Eq. (17al) . one obtains an independent set of ordinary differential 
equations for g^(r,r'-,u), one equation for each £ and u, 

g e (r,r'-,uj) = -<5(r* -r'), (11) 

with 

UM = (i - ™ ) 

Here, r* = r + 2M\n (^7 ~ l) is the radial tortoise coordinate. 

Given two linearly-independent solutions, p(r;cu) and g(r;w), of the ho¬ 
mogeneous version of (fTll) . gi(r,r';uj) is given by 


£(£ + 1) 2 M 

^2 ^ r 3 


( 12 ) 


drl 


+ u 2 - V e (r) 


= P£(r<,o;)g g (r>,a;) 
W(p,q) 


(13) 


5 More generally, Teukolsky EH,!!! showed that the field equations may be separated 
in Kerr spacetime in the gravitational case by making use of the spin-weighted spheroidal 
harmonics in place of the spherical harmonics. A series of works — pioneered by Regge 
and Wheeler [136ll and improved upon by others |l37l - ll43i] — achieved a similar separa¬ 
tion in the Schwarzschild gravitational case by making use of tensor spherical harmonics. 
However, separability alone is not sufficient and there remain some technical issues which 
have yet to be solved before solutions of the Teukolsky or Regge-Wheeler equations could 
be used in a worldline convolution approach. Most important is the issue of gauge; the 
MiSaTaQuWa equations, Eqs. © and were derived in the Lorenz gauge, whereas 
solutions of the Teukolsky and Regge-Wheeler equations are in a gauge different from 
Lorenz gauge. Fortunately, recent work [l 19tl has resolved many of the conceptual issues 
associated with gauge choice. 
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where r> = max(r, r’). r < = min(r, r') and W(p, q ) is the Wronskian. One is 
then faced with the integral over frequencies in the inverse Fourier transform 
appearing in m- This may be achieved through straightforward integration 


Frequency domain Green function, 1=0 mode 




Frequency domain Green function, /=10 mode 




Figure 2: Spherical harmonic modes of the Schwarzschild scalar Green func¬ 
tion. Left: the frequency domain Green function, gi(r,r'',Lo), as a function 
of real frequency for r = r' = 10 M, i = 0 (top) and £ = 10 (bottom). Right: 
the corresponding time domain Green function, gi(r,r'\t — t') computed by 
Fourier transforming the frequency domain Green function using Eq. m 
with cu max = 8.5 (blue, solid line), and by using the time domain methods 
described in Sec. 13.1.31 (orange, dashed line). 


along the real-cu axis (see Fig. ED: the only caveat is that the formally-infinite 
integral over frequencies should be cut off at some finite maximum frequency 
using a smooth window function, for example 


g e (r,r';t- t') = — 


r* 00+26 


9 e(r, r']u)e 


—iu{t—t') /1 _ 


— 00+26 


(1 - erf[2(cu - cj max )])/2du;. 

(14) 

Alternatively, as was done in [69(, the integration contour can be de¬ 
formed into the complex frequency plane following a proposal by Leaver 
144, 1451 ]. In Schwarzschild spacetime, gi{r,r']u) has simple poles in the 
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lower semi-plane at the quasinormal mode frequencies, along with a branch 
cut starting at uj = 0 and continuing along the negative imaginary axis (see 
Fig. [3]). The deformed frequency integral is therefore given by an integral 
over a high-frequency arc, an integral around a branch cut, and an integral 
over the residues at the pole^l. 



Figure 3: Contour deformation in the complex-frequency plane. The 

residue theorem of complex analysis allows one to re-express the integral 
over (just above) the real line of the Fourier modes, ge(r, r / ; uj), as an inte¬ 
gral over a high-frequency arc plus an integral around a branch cut and a 
sum over the residues at the poles. Reproduced from Ref. [69.]. 


The high-frequency arcs can be disregarded as they are likely to only 
contribute at very early times, where the quasilocal approximation can be 
used. There are well-established methods for accurately computing the lo¬ 
cation of the poles (quasinormal mode frequencies) and the residues at the 
poles. The biggest difficulty in the frequency-domain approach is in eval¬ 
uating the branch cut integral, about which very little was known until 
recently (other than asymptotic approximations for e.g. large radius or late 
times). Substantial recent progress has established methods for calculat¬ 


ing this branch cut contribution 14(J ]; these methods were used in 69] to 


compute the self-force in Schwarzschild spacetime. 

The final step in the frequency domain approach is to sum over spherical 
harmonic modes to produce the full Green function. Here, the distributional 
parts of the full Green function can cause poor convergence in the mode- 
sum. Fortunately, there is a straightforward solution to this problem: by 


°In the Kerr spacetime there are indications that there may be additional branch cuts 
to consider 1231. 
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smoothly cutting off the mode sum at large l, 

-i ^max 

G iet (x,x') = — ^2 p e( cos l)ge(r,r';t- t')e ~ £2/n °ut (15) 

rr e=o 

one obtains a mollified retarded Green which is appropriate for use in a 
self-force calculation, and whose sum over £ converges. Empirically, it has 
been found that choosing l cvA ~ ^ max /5 gives good results. 


3.1.3 Time domain approaches 


The frequency domain approach to computing the retarded Green function 
has several shortcomings. It relies on relatively difficult technical calcula¬ 
tions and has poor convergence properties at early times and in scenarios 
where the worldline is not well-represented by a discrete spectrum of fre¬ 
quencies. Recent work has shown that the Green function may also be 
accurately computed (at least for the purposes of self-force calculations) in 
the time domain using straightforward numerical evolutions. A time do¬ 
main calculation sidesteps issues related to a wide frequency spectrum and 
appears to exhibit much better convergence properties at early times. 

There are two closely-related proposals for computing the Green function 


in the time domain. In [7l[], Eq. (I7al) was numerically solved as an initial 
value problem, with the delta-function source approximated by a narrow 
Gaussian. A reformulation of this approach as a homogeneous problem 
with Gaussian initial data was subsequently given in [7(]]. It was found that 
these “numerical Gaussian” approaches are able to approximate the retarded 
Green function remarkably well in a large region of the space required by 
worldline convolutions. The size of the Gaussian limits the scale of the 
smallest features which can be resolved, but it turns out that this is not 
significantly detrimental to a self-force calculation. 

The only regimes where the numerical Gaussian approach is not well 
suited to computing the retarded Green function are at very early and very 
late times. At early times the direct 6(a(x, x)) term (which must be excluded 
from the worldline convolution of the retarded Green function) is smeared 
out and is difficult to isolate from the rest of the Green function. There is no 
conceptual difficulty at late times, but the reality of a numerical evolution is 
that it can only be run for a finite time. Fortunately, both of these issues are 
easily overcome; the former by using the quasilocal expansion at early times, 
and the latter by using a late-tinre expansion of the branch cut integral (see 
Fig. 0]). 
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Schwarzschild retarded Green function along an eccentric geodesic 
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Figure 4: Top: Matched Green function along an eccentric orbit (eccen¬ 
tricity e = 0.5 and semilatus rectum p = 7.2) in Schwarzschild spacetime. 
The full Green function (black) is constructed by combining a quasilocal 
approximation at early times (orange, dashed line) with a time domain cal¬ 
culation at intermediate times (blue, dashed line) and a late-time branch cut 
approximation at late times (green, dashed line). Bottom: Integrating the 
retarded Green function (excluding the direct part) up to some time point 
on the past worldline gives the contribution to the regularized self-field from 
all points on the past worldline up to that time. With an eccentric orbit 
of period ~ 317M, we see that a good approximation to the self-field is 
obtained by including the contributions from less than one full orbit. Note 
that the formal divergences of the Green function are integrable and do not 
cause significant numerical difficulty in computing the self-force. 
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In the time domain approach, each numerical time domain evolution 
gives the Green function G(xo,x') for a single base point xo, and for all 
source points, x. As a result, a single numerical calculation can only be 
used to compute the self-force at a single spacetime point, but it can be 
computed for any past-worldline ending up at that point. This is in contrast 
to other self-force methods, where a single orbit is considered at a time, but a 
single calculation gives the self-force at all points on the orbit. The problem 
of efficiently spanning the parameter spa, ce of base points xo is an ideal 
application of reduced order methods 00. 


3.2 Mode-sum regularization 

The mode-sum regularization scheme, first proposed by Barack and Ori 
.54], has proven highly successful as a computational self-force strategy, 
having been applied to the com put ation of the the self force f or a variety of 
configurations in Schw arzschild 74 -81. 83, 102. 107 - 11G . 113 . 114 . 148 . 149 ] 


and Kerr [89j, |9l|, [92|, 124, ll50l] spacetimes. The success of the method 


hinges on the fact that the first order retarded field diverges in a way which 
can be effectively smoothed out by a spectral decomposition in the angular 
directions. More specifically, the divergence appears as an odd power of 1/s, 
where s 2 = ( g a13 + u a u^)a- a (7-^ is an appropriate measure of distance from 
the worldline. The result is that the 1/s divergence of the held near the 
worldline turns into an infinite sum of modes, each of which are individually 
finite (but possibly discontinuous) on the worldline. The divergence of the 
held then manifests itself through the failure of the (inhnite) sum over modes 
to converge. Conveniently, this odd-in-l/s property of the retarded held also 
holds for its derivatives, both the first derivatives required for the self-force 
and higher derivatives which are useful for computing higher-multipole gauge 
invariants 104 . 1051]. As a result, an arbitrary number of derivatives of the 
hrst order retarded held is represented by an inhnite sum of modes, each of 
which are individually finite (but possibly discontinuous) on the worldline. 
The trade-off is that as the number of derivatives increases the divergence 
of the sum over modes becomes increasingly strong. 

Since the individual modes are finite, numerical calculations of the re¬ 
tarded held and its derivatives can be done in the reduced (t, r ) space with¬ 
out encountering any numerical divergences. The remaining piece of the 
problem is a method for rendering the divergent sum over modes finite. An 
analytic decomposition of the angular dependence of the Detweiler-Whiting 
singular held yields “regularization parameters” which may be subtracted 
mode-by-mode from the numerical retarded held values. Provided all parts 
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of the Detweiler-Whiting singular field which don’t vanish on the worldline 
are included in the calculation, the regularized sum over modes is convergent 
on the worldline and one never encounters any numerical divergences. 

It is important to note that the use of the Detweiler-Whiting singular 
field is not merely a convenience; it also provides well-motivated physical 
grounds to justify the mode-sum regularization approach. Other ad-hoc ap¬ 
proaches based on identifying the asymptotically divergent contributions to 
the mode-sum often produce the correct regularization parameters, but their 
use is difficult to justify rigorously and can easily lead to incorrect results. 
Crucially, there is no way of knowing whether such an ad-hoc regularization 
procedure is producing a physically correct result, or if important contribu¬ 
tions are being overlooked. They should therefore not be relied on without 
extreme care. 

Unfortunately, despite its success in first order calculations, mode-sum 
regularization alone is not an effective tool for computing the second order 
self-force. The reason for this is straightforward: the second order retarded 
field contains divergent terms which appears in the form of even powers of 
1/s. Intuitively, this arises from the fact that the second order field contains 
terms which depend quadratically on the first order field. Unlike the odd- 
in-l/s case, the angular decomposition of 1/s 2 leads to individual modes 
which diverge logarithmically as the worldline is approached. Fortunately, 
all is not lost for the mode-sum method as a computational tool at second 
order. Provided the leading order logarithmic divergence is subtracted by 
some other means (for example, using the effective source approach), reg¬ 
ularization parameters may be used to accelerate the rate of convergence 
of the mode sum. Such a hybrid scheme complements the generality of the 
effective source approach with the computational efficiency of mode-sum 
regularization. 

3.2.1 Regularization parameters 

The computation of regularization parameters has been addressed by a series 
of calculations stemming from the original Barack-Ori derivation. Barack 
and Ori’s original work gave the first two self-force regularization param¬ 
eters for scalar, electromagnetic and Lorenz-gauge gravitational cases in 
both Schwarzschild and Kerr spacetimes. These are sufficient for computing 
the regularized self-force, but they yield mode-sums which have relatively 
poor quadratic convergence with the number of modes included. Since the 
Detweiler-Whiting regularized field is a smooth, homogeneous function in 
the vicinity of the worldline, its mode-sum representation converges expo- 
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nentially. This spectral convergence is spoiled by the fact that only the 
portion of the Detweiler-Whiting singular field which does not vanish on 
the worldline is subtracted by the leading order regularization parameters. 
The mode decomposition effectively contains information about the exten¬ 
sion of the field off the worldline to the entire two-sphere, so the regularized 
field contains residual pieces of the Detweiler-Whiting singular field off the 
worldline, but on a two-sphere of the same radius. By deriving higher-order 
regularization parameters one can subtract these residual pieces order-by¬ 
order, leaving a mode-sum which is more and more rapidly convergent (see 
Fig- E]). The derivation of these higher-order parameters is closely related to 
the computation of the quasilocal expansion of the Green function and has 
been addressed in a series of papers: the first higher order parameter was 
given in 76] for the case of a scalar charge on a circular orbit in Schwarzschild 
spacetime, and for eccentric geodesic orbits in [8o| . This was subsequently 
extended by several further orders for equatorial geodesic motion in the 
scalar, electromagnetic and gravitational cases in both Schwarzschild [§3] 
and Kerr 90] spacetimes. Recent work has also produced parameters for 


accelerated worldlines [72, 93]. 



Figure 5: Spherical harmonic modes of the self-force for a point scalar 
charge on a circular orbit of radius ro = 6 M in Schwarzschild spacetime. By 
subtracting analytically determined high-order “regularization parameters”, 
the sum over modes is rendered more and more rapidly convergent. Each 
regularization parameter, F^\ behaves asymptotically for large-£ as 1 /i n . 
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Within these derivations of high-order regularization parameters, a sub¬ 
tle ambiguity appears through the elevation of the four-velocity from a quan¬ 
tity defined on a worldline to a quantity defined everywhere on the two- 
sphere. A natural, covariant choice is to define this off-worldline extension 
through parallel transport, u a = g a a u a . However, in practical calcula¬ 
tions it is often convenient to make a coordinate choice. For example, a 
common choice is to define the extension in terms of “constant coordinate 
components”, i.e. to define u a such that its components in some coordinate 
system have a constant value everywhere on the two-sphere. This is a per¬ 
fectly valid choice, as is any other choice where u a agrees with the actual 
four-velocity when evaluated on the worldline. The only caveat is that the 
regularization parameters beyond the leading two orders change depending 
on the particular choice of extension. As a result, in order to use higher- 
order regularization parameters it is essential that a compatible choice of 
off-worldline extension is used in the retarded field calculation. 


3.2.2 Choice of basis and gauge 


There are two other factors which must be considered in the mode-sum 
scheme: the choice of angular basis functions for the spectral decomposition 
and the choice of gauge in the electromagnetic and gravitational cases. 

The choice of basis functions is typically motivated by their ability to 
produce separability of the retarded field equations; an appropriate choice of 
basis functions yields an independent set of equations in (t, r) space for each 
individual mode. For the case of a scalar charge in Schwarzschild space- 
time, the appropriate basis functions are the standard spherical harmon¬ 
ics. In the Kerr case the spheroidal harmonics are required for separability. 
For electromagnetic and gravitational cases, there are several choices. For 
Schwarzschild spacetime one can choose between a tensor harmonic basis and 
a basis of spin-weighted spherical harmonics. In the Kerr case o nly t he spin- 


weighted spheroidal harmonics are known to yield separability 134 . 135 ]. 


This choice of angular basis also affects the regularization parameters. 
The parameters for scalar spherical harmonic modes, tensor harmonic modes, 
spin-weighted spherical harmonic modes and spheroidal harmonic modes are 
all potentially different. However, it is always possible project the tensor, 
spin-weighted, or spheroidal harmonic modes onto the scalar spherical har¬ 
monic basis. In fact, since the regularization parameters are most easily 
obtained for the scalar spherical harmonic basis, calculations of the retarded 
field are typically done using the computationally most convenient basis and 
the result is then projected onto the scalar spherical harmonics before the 
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regularization and sum-over-modes steps are done tm. 


A more difficult issue in the electromagnetic and gravitational cases is the 
choice of gauge. The Detweiler-Whiting singular field is defined in Lorenz 
gauge (since it is derived from a Lorenz gauge Green function), but numer¬ 
ical calculations of the retarded field are more easily done in either radia¬ 
tion or Regge-Wheeler gauge. The existence of tensor spherical harmonics 
makes a Lorenz-gauge calculation possible, if somewhat cumbersome in the 


Schwarzschild case 111, 118 . 120 l |. One obtains a coupled set of 10 equations 


for the metric perturbation, but there is no coupling between different differ¬ 
ent tensor-harmonic modes. Unfortunately, this does not hold for the Kerr 
case as there are no known tensor spheroidal harmonics. Rather than trying 
to derive tensor spheroidal harmonics, a better approach is to work with the 
relatively straightforward Teukolsky equation in radiation gauge 124]. The 
difficulty then is in identifying the appropriate regularization parameters, 
particularly since the gauge transformation from Lorenz gauge is itself often 
singular. This remained an open problem until recently, when an under¬ 
standing of how to apply the mode-sum regularization scheme in radiation 
gauge was finally established jll9] and applied in a self-force calculation 


108|]. 


3.2.3 Mode-sum regularization in the frequency domain 

The mode-sum scheme provides a method for computing the regularized 
self-force using numerical data for the modes of the retarded field in the 
reduced ( t , r) space. There are, however, several possibilities for computing 
this numerical data. One option relies on a further decomposition of the 
time dependence into Fourier-frequency modes, in an analogous way to the 
frequency-domain Green function described in Sec. 13.1.21 above. Using a 
spin-weighted spheroidal harmonic basis, this leads to a radial equation for 
the Teukolsky function, 


A~ s 4- (a s+1 ^) + V(r)R = T(r) 
dr V dr / 

with the potential given by 


V(r) 


K 2 - 2 is(r - M)K 
-^---1- 4 irus - A. 


(16) 


(17) 
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Here, A = r 2 — 2 Mr + a 2 , K = (r 2 + a 2 )w — am, A is an eigenvalue of the 
spheroidal equation, 


1 d 

sin 9 dO 


( sin 9-^-\ — a 2 cu 2 sin 2 6 — 
V dd J 


(m + s cos 9)^ 


sin 


— 2acus cos 0 + s + 2maw + A] s S^ m = 0 (18) 


and a and M are the Kerr spin and mass parameters. For a = 0, A = 
(£ — s)(£ + s + 1) and this reduces to the Schwarzschild wave equation de¬ 
composed into spin-weighted spherical harmonics (in radiation gauge for the 
gravitational case), while for s = 0 it reduces to the scalar wave equation. 
A decomposition of the point-particle source term yields a source involving 
5(r — r p ) (and in some cases its derivative), where r p is the radial loca¬ 
tion of the particle. In practice, the frequency domain mode-sum approach 
proceeds in the same way as for the Green function: two independent solu¬ 
tions of the homogeneous radial equation are obtained and are matched at 
the particle’s locatior@. Then, the distributional sources do not introduce 
any numerical difficulty as they simply appear as jumps when matching the 
homogeneous inner and outer solutions. 

The solutions of the Teukolsky equation for a given {£, m, u) can be 
obtained either through straightforward numerical integration of the radial 
ordinary differential equa tion (typically with some modifications to improve 
numerical accuracy 
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152]) or as an approximation in the form of an 
infinite convergent series of hypergeometric functions. The latter method 


is based an idea originally developed by Leaver 1531 ] and now commonly 


referred to as the “MST” method, after Mano, Suzuki, and Takasugi 154] 


who reformulated it into its current form. It provides an efficient and highly 
accurate method for computing solutions of the radial equation. For ex¬ 
ample, recent results have used it to compute solutions accurate to several 
hundred decimal places, allowing the solutions to be used to determine pre¬ 


viously unknown high-order post-Newtonian parameters 155— 1591] . For a 


comprehensive review of the MST approach, see the living review by Sasaki 
and Tagoshi Q and references therein. 

The frequency-domain approach is particularly appropriate in scenarios 
where the worldline is well represented by a discrete spectrum involving a 
small number of frequencies. In such cases, mode-sum regularization is by 
far the optimal choice and is unparalleled in its accuracy and computational 


'For eccentric orbits where the particle can not be considered to be at a single ra¬ 
dial location in the frequency domain this matching must be modi fied slightly using, for 
example, the method of extended homogeneous solutions ms, mi- 
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efficiency 16l|. The prototypical example is a circular orbit, in which only 
a single frequency is present. In the case of eccentric equatorial orbits (and 
inclined circular orbits in the Kerr case), there are two fundamental fre¬ 
quencies, and also an infinite number of higher harmonics produced from 
combinations of the fundamental frequencies. For mildly eccentric orbits 
this does not cause a great deal of difficulty. However, for more eccentric or¬ 
bits (with eccentricities e > 0.5) an increasingly large number of frequencies 
must be included and the competitive advantage of the frequency domain 
approach is lost 161]. Even worse, generic geodesic orbits in Kerr space¬ 


time have three fundamental frequencies and the computational difficulty is 
so high that a calculation has yet to be attempted. Similarly, unbound orbits 
and other cases such as radial infall are not well suited to frequency domain 
methods. Apart from these deficiencies, the frequency domain mode-sum 
approach has been highly successful for producing results. 


3.2.4 Mode-sum regularization in the time domain 

The mode-sum scheme may also be applied in the time domain by skipping 
the Fourier decomposition step and instead solving a set of 1 + ID partial 
differential equations in (t, r ) space. The main difficulty then is in appropri¬ 
ately handling the distributional s ourc e term which has the form S(r — r p (t)). 


113. lUlhfii]. is to use a discretised rep- 


One solution, used in [4] 
resentation of a delta function and to construct the computational grid such 
that the worldline only ever passes between grid points. 

An alternative approach is to reformulate the problem in an analo¬ 
gous way to the frequency-domain method. By splitting the computa¬ 
tional domain into two domains — one either side of the particle — the 
delta-function source can be reformulated in terms of a jump in the fields 
and their derivatives at the interface of the two domains. This method is 
well suited to highly-accurate spectral numerical methods as all of the nu¬ 
merically evolved fields are smooth fun ction s. It was implemented using 
discontinuous-Galerkin methods in lod . 115 ] and using Chebyshev pseudo- 


spectral methods in |78l . 81]. Eccentric orbits present a small additional 


complexity in this case as the particle must always lie on the domain in¬ 
terface. This is easily achieved by introducing a time-dependent mapping 


between the computational and physical coordinates of the system 81, 1061]. 


23 



























3.2.5 Limitations of the mode-sum regularization scheme 

Despite its resounding success to date, the mode-sum regularization scheme 
has some unfortunate disadvantages which make it ill-suited as a general- 
purpose method for self-force calculations: 

1. Its application in the Kerr case is only straightforward in the frequency 
domain, since the field equations are not separable in the time domaiiu 

2. It relies on the use of Lorenz gauge for regularization in the gravita¬ 
tional case. This is not a major issue in the Schwarzschild case since 
the tensor spherical harmonics may be used to decouple the Lorenz 
gauge field equations in the angular directions (leaving 10 coupled 
1 + ID equations for each £,m mode). However, there are no know 
tensor spheroidal harmonics which would be required for the Kerr case, 
and even if there were they would likely not be applicable in the time 
domain (again, it is conceivable that a coupled system of equations 
involving tensor spherical harmonics could be used, but the coupling 
would result in considerable complexity). 

3. It is not applicable beyond first perturbative order, since the modes of 
the second order perturbation diverge logarithmically near the world- 
line. 


The first two issues are not necessarily showstoppers. There have been 
several attempts at 1 + ID time-domain implementations using coupled 


spherical-harmonic modes in the Kerr case 163. Il64f ]. and recent progress 


on reformulating mode-sum regularization for radiation gauge has clarified 
the regularization issue [119|. The third point, however, appears insur¬ 


mountable. For these reasons, among others, the effective source method, 
described in the next section, was developed. 


3.3 Effective source approach 

Proposed in 2007 as a solution to the shortcomings of mode-sum regulariza¬ 
tion, the effective source approach provides an alternative method for han¬ 
dling the divergence of the retarded field. Rather than first computing the 

8 It may still be possible to use mode-sum regularization for the Kerr case in the time 
domain by decomposing the field equations into spherical harmonics and evolving the 
resulting infinitely coupled set of 1 + ID partial differential equations in a similar manner 
to the Schwarzschild case. 

9 Note that the effective source proposed by Lousto and Nakano [65| is similar in spirit, 
but differs in that it is not derived from the Detweiler-Whiting singular field. 
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retarded field and then subtracting the singular piece as a post-processing 
step, one can instead work directly with an equation for the regular field. 
This idea — independently proposed by Barack and Golbourn [56j and by 
Vega and Detweiler [571] — has the distinct advantage of involving only reg¬ 
ular quantities from the outset, making it applicable in a wider variety of 
scenarios than the mode-sum scheme. In particular, since it does not rely 
on a mode decomposition of the retarded field, it can be used by any nu¬ 
merical prescription for solving the retarded field equations, whether in the 
frequency domain (where the method is really just a generalisation of mode- 
sum regularization) or in the time domain as a 1 + ID, 2 + ID or even 3 + ID 
problem. 

The basic idea is to use the split of the retarded field into regular and 
singular pieces, Eq. ([4|l . to rewrite the field equations, Eqs. (flail . (|2aj) and 
fall as equations for the Detweiler-Whiting regular field, 


(□ - £i?)$ R = (□ - eR)($ ret - 4> s ) 

= -4irp - (□ - £i?)4> s , 


(19a) 


(□5/ - Ri)A f = (06J - Ri)(Af - A%) 

= -47re J g aa iu a '8 ±(x,z(t')) dr' - (05^ - (19b) 

(□<WV + 2C-aV)^ = (D<5aV + {?)(}£$ ~ h%) 

=- 167 TH J g a ^ a u a '5 A (x,z(T'))d,T' - (U5c?5p 5 + 20^ p S )(h^ s ). 

(19c) 

If the singular held used in the subtraction is exactly the Detweiler-Whiting 
singular held, then the two terms on the right hand side of this equation 
cancel and the regularized held would be a homogeneous solution of the 
wave equation. Unfortunately, one typically does not have access to an exact 
expression for the singular held. Indeed, the Detweiler-Whiting singular held 
is only defined through a Hadamard parametrix which is not even defined 
globally. Instead, the best one can typically do is a local expansion which 
is valid only in the vicinity of the worldline. Borrowing the language of 
Barack and Golbourn, we refer to an approximation to the singular held as 
a “puncture” held, <f> s ss <h p , ^4® ~ h^g ~ Then, the corresponding 
approximate regular held — referred to as the “residual held” — is no longer 
a solution of the homogeneous equation, but instead is a solution of the 
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sourced equation with an effective source (see Fig. [6]) which is defined to be 
the right-hand side of Eq. ()19D with the puncture field substituted for the 
singular field 


S eS = —Airq J 5 4 (x, z(r'))dT' - (□ - £R)<S> P , (20a) 

= -47re J g aa >u a 'S 4 (x, z{t'))(It' - {U5 a p - R^)(A^), (20b) 

Sfp = -16vr h J g a ,( a u a 'gp^u 13 '8 4 (x, z(t'))(It' - (D^ 7 ^ 5 + 2 C a ' y p S )(h p s ). 

(20c) 

Note that the presence of a distributional component of the source on 



Figure 6: The effective source for a scalar particle on a circular orbit 
in Schwarzschild spacetime. The source is generically non-smooth in the 
vicinity of the worldline (left), but the smoothness can be improved by 
incorporating higher-order parts of the Detweiler-Whiting singular field into 
the source. 

the worldline is merely a formal prescription; in practice the puncture field 
is chosen so that it exactly cancels this distributional component on the 
worldline. This effective source is then finite everywhere, but has limited 
differentiability on the worldline. This makes it well-suited to numerical 
implementations since no divergent quantities are ever encountered. The 
only numerical difficulty arises from the non-snroothness of the source in 
the vicinity of the worldline (see Fig. [6]). which leads to numerical noise in 
the computed self force. The noise can be reduced by making the source 
smoother using higher-order parts of the Detweiler-Whiting singular field. 
As shown in Fig. 0 at the same numerical resolution a higher order source 
(C 2 in this case) eliminates the vast majority of the numerical noise that 
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is present when using a lower order source (C° in the case in the figure). 
The cost of this improved accuracy is a source which is considerably more 
complicated, and costly to compute in a numerical code. 


Scalar self-force for e=0.5, p=7.2 geodesic 



Figure 7: Radial self-force for a scalar particle on an eccentric orbit in 
Schwarzschild spacetime, computed with a 3 + ID implementation of the 
effective source scheme. Here, the independent variable x is a “relativis¬ 
tic anomaly parameter” defined through r = pM/( 1 + ecosx)> with p the 
semilatus rectum and e the eccentricity. The high-frequency errors using a 
continuous but non-differentiable source (blue) are dramatically decreased 
by using a twice differentiable source (orange) obtained from a higher-order 
approximation to the Detweiler-Whiting singular field. For reference, a 
highly accurate value computed using frequency-domain mode-sum regular¬ 
ization is also included (dashed, black). Figure based on version presented 


in Ref. 165l ]. 


An additional level of complexity arises from the fact that the puncture 
field is defined only in the vicinity of the worldline. To avoid ambiguities 
in its definition far from the worldline, one must ensure that the puncture 
field goes to zero there. This is most easily achieved by multiplying it by a 
window function, W, with properties such that it only modifies the puncture 
field in a way that its local expansion about the worldline is preserved to 
some chosen order. In a first-order calculation of the self-force, it suffices to 
choose W such that W(x p ) = 1, W'{x p ) = 0, W"{x p ) = 0 and W = 0 far 
away from the worldline. The residual field (see Fig. [8]) then obeys 
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Figure 8: The i = l,m = 1 spherical-harmonic mode of the residual 
field, <h res , for a scalar point particle on a circular orbit in Schwarzschild 
spacetime. This was produced in J86(] using a frequency-domain effective 
source approach. 


(□ - £i?)$ res = S eS (21a) 

(asj - R? a )A) T = Sf (21b) 

(□£a v + icjffyqi = sfp ( 21 c) 

and has the properties 

$ res (x p ) = <I> R (x p ), V Q <f> res (x p ) = V a <f> R (x p ), 

T res (x) = $ ret (x) for x £ supp(W), (22a) 

A™(x p ) = A R (x p ), V a A™(x p ) = V a A R (x p ), 

A™ s (x) =A™\x) for x ^ supp(W), (22b) 

K$( X p) = h aP ( x p), V Q /l^(x p ) = V a h^( Xp ), 

hapix) = haH x ) for x £ supp(W). (22c) 


As the residual held coincides with the retarded held far from the par¬ 
ticle we can use the usual retarded held boundary conditions when solving 
Eqs. (I21al) . (I21bl) and (I21cl) . The details of a numerical implementation of 
the effective source approach are then much the same as for the Green func¬ 
tion and mode-sum regularization schemes. One can use either a frequency 
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domain or a time domain method for solving the field equations, the key 
differences now being that there is no restriction to 1 + ID, and that the 
effective source must be included as a source term. A more thorough re¬ 
view of the effective source approach — including a detailed description of 
methods for ob tain ing explicit expressions for the effective source — can be 


found in Refs. 126|, Il66( ] 


The effective source approach has been successfully applied in the fre¬ 
quency domain 18611. and in the time domain in 1+1D [57]], 2+1D |8^,i96j, 122] 
and 3 + ID [85j, 87, 88] contexts. It has also been formulated — but not 
yet implemented — in the gravitational case to second order in the mass 
ratio {3, 0; for the second order problem it is currently the only viable 
computational strategy. Since the effects of the second-order metric per¬ 
turbation will be very small — being suppressed by two orders of the mass 
ratio relative to test body effects — it is likely that a highly accurate nu¬ 
merical scheme will be required, suggesting a frequency domain treatment 
of the problem where one encounters ordinary differential equations (ODEs) 
which are relatively easy to solve numerically to high accuracy. 


4 Evolution schemes 

The calculation of the self-force is only the first stage in the production of 
an inspiral model. Another critical component is a scheme for evolving the 
orbit using this self-force information. Various approaches have been pro¬ 
posed, each of which brings with it its own advantages and disadvantages. 
Approaches which make approximations in the self-force used to drive the 
inspiral can give substantial decreases in the computational cost of an inspi¬ 
ral simulation, but come at the cost of ignoring potentially relevant physical 
effects. 

4.1 Dissipation driven inspirals 

A straightforward model for the inspiral can be obtained from energy bal¬ 
ance considerations. Using a flux calculation — in which fluxes of energy and 
angular momentum are obtained by evaluating the point-particle retarded 
field near the horizon and out at infinity - the entire issue of regularization is 
avoided and one obtains an approximation to the contributions to the inspi¬ 
ral coming from dissipative self-force effects. However, this is inadequate for 
capturing all relevant effects from the first-order self-force. Being a dissipa¬ 
tive approximation, it completely misses all conservative corrections to the 
motion. Furthermore, there is no well-defined way of associating fluxes far 
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from the worldline with an instantaneous local self-force on the worldline; a 
flux model can only be used to drive an inspiral in a time-averaged sense, 
ignoring some potentially important dissipative contributions. Nevertheless, 
the relative straightforwardness and computational efficiency of their imple¬ 
mentation have made flux models a compelling approach to assessing qual¬ 
itative features of self-force driven inspirals. These “kludge” models have 
been used to produce kludge waveforms for EMRI systems which capture at 
least some of the qualitative physical effects 44, |45|, |48|, [491. 11621 ]. 




Figure 9: Radial (left) and time (right) components of the self-force through 
one radial cycle, for three different geodesic orbits in Schwarzschild space- 
time. Solid lines indicate the full self-force and dashed lines indicate the 
conservative-only (left) or dissipative-only (right) pieces. Arrows denote the 


direction of the geodesic motion. Reproduced from Ref. 1871 ]. 


To improve on the flux-based dissipative model, one can instead make use 
of a half-retarded minus half-advanced scheme 53 , 13-0 to compute a local, 
instantaneous dissipative self-force (see Fig. [9]). This approach captures all 
dissipative effects responsible for driving the inspiral, but still neglects small 
but potentially important conservative corrections to the orbital phase of the 
system 


121 |] . 


4.2 Osculating self-forced geodesics 

In order to account for both conservative and dissipative effects from the 
first order self-force, it is essential to build an inspiral model around a full 
calculation of the local self-force. This, however, still leaves some flexibility 
in the choice of method for computing this local self-force. One option, based 
on the osculating geodesics framework [1671 ] is particularly compelling as it 
allows for dramatic improvements in computational efficiency by separat- 
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ing the coupled problem of simultaneously solving for the orbit and for the 
regularized retarded field into two independent, largely uncoupled compu¬ 
tational problems. The basic idea is to perturbatively expand the worldline 
about a geodesic of the background spacetime, z(r) = ~o(t) + ‘ • • . Then, the 
self-force at first order is a functional not of the evolving orbit z(t), but of 
the geodesic orbit, zo (r), which is instantaneously tangent to the worldline. 
This expansion is valid in the adiabatic regime, where the orbit is evolving 
slowly and the difference between the geodesic and evolving orbits is small; 
the error introduced by the approximation appears at the same order in the 
equations of motion as a second-order perturbative correction to the field 
equations. 

The improvements in computational efficiency brought about by the 
osculating geodesics approximation are dramatic. The problem of self- 
consistently computing the regularized self-force coupled to an arbitrarily 
evolving orbit is reduced to the much simpler problem of determining the 
self-force for geodesic worldlines. Since orbits of black hole spacetimes are 
parametrised by at most three conserved quantities (energy and angular 
momentum in the Schwarzschild case, supplemented by the Carter constant 
in the Kerr case), it is computationally tractable to span the entire pa¬ 
rameter space of moderately-eccentric geodesic orbits. Even better, the 
highly-accurate frequency domain mode-sum method can be used because 
bound geodesic orbits are efficiently represented by a small frequency spec¬ 


trum. In Ref. 1211 ] this approach was explored over an entire radiation- 


reaction timescale for moderately-eccentric inspirals in the Schwarzschild 
spacetimc0. This was achieved by tabulating the values for the self-force in 
a relevant portion of the energy-angular momentum phase space of geodesic 
orbits, and using an interpolated model of the tabulated results to drive an 
orbital inspiral. 


4.3 Self-consistent evolution 

Unfortunately, the adiabatic approximation responsible for the dramatic im¬ 
provements in computational efficiency brought by the osculating geodesics 
framework also introduces errors in the equation of motion at second pertur¬ 
bative order, making the method inadequate for the purposes of precision 
EMRI astrophysics. One possible solution, implemented in 88] for the scalar 
field case, is to avoid expanding out the worldline and instead evolve the 
self-consistent coupled system of equations, Eqs. (fTa|) and (llbl) : (l2a|) . ()2bl) 


10 See also Ref. 16 
circular orbits. 


for an approximate version which assumes a sequence of quasi- 
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and ([2c]) : or (l3al) and (I3b|) . This is a computationally much more difficult 
problem as the field equations must be solved in the time domain and one 
cannot rely on an efficient off-line tabulation of self-force values. This self- 
consistent evolution can in principal be implemented using a time-domain 
mode-sum scheme, but in practice implementations have instead used the 
effective source approach since that will be necessary at second perturbative 
order (see Sec. 13.21 for an explanation of this issue). 

While a self-consistent evolution incorporates all effects contributing to 
the first order self-force, and only neglects second order contributions from 
the second order field equations, this comes at the cost of computational effi¬ 
ciency. Whereas an osculating geodesics framework can evolve a large num¬ 
ber (~ 10,000) of orbits with ease once the initial off-line tabulation phase is 
complete a the self-consistent scheme requires a long calculation of the 
solutions of the first-order field equations for each new orbit. In reality, this 
limits the practicality of the scheme to tens of orbits using current methods. 
The self-consistent evolution scheme is therefore most valuable as a bench¬ 
mark against which other, less accurate but more efficient methods can be 
validated. In fact, comparisons between the self-consistent and osculating 
geodesics scheme for the scalar case indicate that the osculating geodesics 
scheme performs remarkably well [ 1691 ]. 


4.4 Two-timescale expansions 

The osculating geodesics scheme is fast, but inaccurate, while the self- 
consistent evolution scheme is accurate, but slow. This begs the question of 
whether there is a middle ground which incorporates most of the accuracy 
of a self-consistent evolution while maintaining the computational efficiency 
of the osculating geodesics scheme. One promising possibility is that the 
use of a two-timescale expansion could give just such a scheme. In the 
two-timescale scheme, rather than using an expansion about a background 
geodesic (which is valid over short timescales characterised by the orbital 
motion), one introduces an additional radiation reaction timescale into the 
problem and incorporates the relevant effects over this radiation reaction 
timescale. 

The relevant two-timescale expansion of the worldline equations of mo¬ 
tion was completed in Ref. [0], and follow-up work has improved our under¬ 


standing of important resonance effects for inspiral orbits 0, 00. In 


order to consistently incorporate a two-tinrescale expansion into an orbital 
evolution scheme, it will also be necessary to have a two-time expansion 
of the field equations. Promising progress towards such an expansion was 
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recently reported in 18 o( ]. indicating that it is likely that a self-consistent 
orbital evolution using a two-timescale expansion is indeed feasible. 


5 Discussion 


The techniques described in this review represent the current state of the 
art of self-force calculations. The three primary approaches: worldline con¬ 
volutions, mode-sum regularization and the effective source approach can be 
considered complimentary, with each having regimes where they are most 
appropriate: 

• The mode-sum scheme gives unparalleled accuracy (particularly in the 
frequency domain) for orbits which have a small frequency spectrum. 

• The worldline convolution method provides valuable physical insight 
and can be easily applied to arbitrary orbital configurations, including 
those inaccessible by other means. 

• The effective source approach can be used in arbitrary spacetimes with¬ 
out relying on any symmetries, and also stands out as the most appli¬ 
cable to a second order calculation. 


There still remain important developments to be made in each case. For 
example, the mode sum and effective source schemes are still under develop¬ 
ment for the Kerr gravitational case, and worldline convolution approaches 
have yet to be fully applied to the gravitational case for any spacetime. 

Finally, it should be pointed out that this review is not an exhaustive 
exposition of all self-force computation strategies. Many other calculations 
have not been described, including alternative regu larization strategies [66 


181 


182l |. near-horizon waveform calculations 1831 - 186 ], methods based on 


effective f ield theo ry 187 - 19 til ] . analytic calculations in particularly simple 


cases 


95, 191, 192], black holes in higher dimensions 


194 


193], and calculations 


195 ] and cosmological 


in non- black hole spacetimes such as wormholes 

models 19(1 197 ], More details can be found in the reviews 60|, l6l|, and 
references therein. 
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